/** Fichier implémentant une classe de production de bruit cohérent.
 *
 * Ce fichier contient l'implémentation de la classe Bruit.
 *
 * @file Bruit.cpp
 * @author Johan "Solo" Luisier
 *
 * @date 2010
 */
#include "Bruit.hpp"

#include <cmath>

#include <boost/assign/list_of.hpp>

#include "Utils/ATL.hpp"
#include "Utils/borne.hpp"

using namespace std;
using namespace Perlin;
using namespace boost::assign;

const double Bruit::normalisation = pow( 2., sizeof( int ) * 8 ) - 1.;

const double Bruit::Pi = acos( -1. );

/**
 * Constructeur standard. Initialise les différents champs de la classe.
 *
 * @param periode période de la perturbation fondamentale.
 * @param pers rapport entre les amplitudes de deux harmoniques.
 * @param oct nombre d'harmoniques à utiliser.
 * @param fct type d'interpolation à utiliser.
 */
Bruit::Bruit(const unsigned& periode, const double& pers,
	     const unsigned& oct, const InterpolFct& fct)
  : Fonction( fct ), Periode( periode ), dPeriode( periode ),
    Persistance( pers ), Octaves( oct )
{}

/**
 * Le destructeur n'a rien de spécial à faire.
 */
Bruit::~Bruit()
{}

/**
 * Permet de changer le type d'interpolation à utiliser.
 *
 * @param fct nouveau type d'interpolation à utiliser.
 */
void Bruit::fonction(const InterpolFct& fct)
{
  Fonction = fct;
}

/**
 * Permet de changer la période de la perturbation fondamentale.
 *
 * @param t nouvelle valeur de la période de la perturbation
 * fondamentale.
 */
void Bruit::periode(const int& t)
{
  Periode = t;
}

/**
 * Permet de changer le rapport entre les amplitudes de deux
 * harmoniques successives.
 *
 * @param pers rapport entre les amplitudes de deux harmoniques.
 */
void Bruit::persistance(const double& pers)
{
  Persistance = pers;
}

/**
 * Permet de changer le nombre d'harmoniques à utiliser lors de la
 * génération du bruit.
 *
 * @param oct nouvelle valeur du nombre d'harmoniques à utiliser.
 */
void Bruit::octaves(const unsigned& oct)
{
  Octaves = oct;
}

/**
 * Permet de générer la valeur du bruit pour des entiers en une, deux
 * ou trois dimensions.
 *
 * La dimensionalité est directement déterminée d'après le nombre
 * d'arguments.
 *
 * @param x graine 1.
 * @param y graine 2.
 * @param z graine 3.
 */
double Bruit::operator()(const int& x, const int& y, const int& z) const
{
  double resultat(0.), norme( 0. );
  if ( y == numeric_limits< int >::max() )
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ) ) *
	    pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
  else if ( z == numeric_limits< int >::max() )
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ),
				      y * pow( 2, i ) ) *
	    pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
  else
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ),
				      y * pow( 2, i ),
				      z * pow( 2, i ) ) *
	    pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
   
  return borne( resultat / norme, 0., 1. );
}

/**
 * Permet de générer la valeur du bruit pour des nombres réels en une,
 * deux ou trois dimensions.
 *
 * La dimensionalité est directement déterminée d'après le nombre
 * d'arguments.
 *
 * @param x graine 1.
 * @param y graine 2.
 * @param z graine 3.
 */
double Bruit::operator()(const double& x, const double& y,
			 const double& z) const
{
  double resultat(0.), norme( 0. );
  if ( y == numeric_limits< double >::max() )
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ) )
	    * pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
  else if ( z == numeric_limits< double >::max() )
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ),
				      y * pow( 2, i ) )
	    * pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
  else
    {
      for (unsigned i(0); i < Octaves; i++)
	{
	  resultat += bruitInterpole( x * pow( 2, i ),
				      y * pow( 2, i ),
				      z * pow( 2, i ) )
	    * pow( Persistance, i );
	  norme += pow( Persistance, i );
	}
    }
 
  return borne( resultat / norme, 0., 1. );
}

/**
 * Génère un nombre dans l'intervalle [0;1] à partir de l'entier donné
 * en argument.
 *
 * La fonction génère un nombre entier par permutation, qui est
 * ensuite transformé en un nombre réel dans [0;1].
 *
 * @param x "graine".
 *
 * @return valeur de la valeur de la permutation normalisée au plus
 * grand entier.
 */
double Bruit::bruit(const int& x) const
{
  int var( ( x << 13 ) ^ x );
  unsigned resultat( var * var * 15985 + 735971 );
  resultat *= var;
  resultat += 675332129;
  resultat ^= 0x7FFFAAFF;
  double ret( resultat );
  ret /= normalisation;
  return ret;
}

/**
 * Permet de générer un nombre dans l'intervalle [0;1] à partir des
 * deux nombre entiers donnés en arguments.
 *
 * La fonction appelle en fait le générateur de bruit en 1 dimension.
 *
 * @param x "graine" 1.
 * @param y "graine" 2.
 * @return valeur de la valeur de la permutation normalisée au plus
 * grand entier.
 */
double Bruit::bruit(const int& x, const int& y) const
{
  return bruit( static_cast<int>( bruit( y ) * normalisation ) *
		( x + 1 ) ) ;
}

/**
 * Permet de générer un nombre dans l'intervalle [0;1] à partir des
 * trois nombre entiers donnés en arguments.
 *
 * La fonction appelle en fait le générateur de bruit en 2 dimensions
 * et en 1 dimension.
 *
 * @param x "graine" 1.
 * @param y "graine" 2.
 * @param z "graine" 3.
 * @return valeur de la valeur de la permutation normalisée au plus
 * grand entier.
 */
double Bruit::bruit(const int& x, const int& y, const int& z) const
{
  return bruit( static_cast<int>( bruit( x, y ) * normalisation ) * ( z + 1 ) );
}

/**
 * Méthode permettant d'interpoler la valeur du bruit en une
 * dimension. Le bruit est donné par plusieurs éléments de la
 * permutation (générés par bruit).
 *
 * Pour une interpolation linéraire ou cosinus, il faut deux points
 * (l'entier précédent et l'entier suivant), pour une interpolation
 * cubique, il faut utiliser quatre points.
 *
 * @param pts valeur des bruits pour les entiers entourant la valeur
 * pour laquelle on veut la valeur du bruit.
 * @param x valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 *
 * @return la valeur interpolée du bruit.
 */
double Bruit::inter(const vector< double >& pts, const double& x) const
{
  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      if ( pts.size() != 2 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	return interLin1D( pts, x );
    case Cubique:
      if ( pts.size() != 4 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	return interCub1D( pts, x );
    default:
      cerr << "Fonction inconnue '" << Fonction << "'" << endl;
      return 0.;
    }
}

/**
 * Méthode permettant d'interpoler la valeur du bruit en deux
 * dimensions. Le bruit est donné par plusieurs éléments de la
 * permutation (générés par bruit).
 *
 * Pour une interpolation linéraire ou cosinus, il faut deux points
 * (l'entier précédent et l'entier suivant), pour une interpolation
 * cubique, il faut utiliser quatre points.
 *
 * @param pts valeur des bruits pour les entiers entourant la valeur
 * pour laquelle on veut la valeur du bruit.
 * @param x valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 * @param y valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 * @return la valeur interpolée du bruit.
 */
double Bruit::inter(const vector< vector< double > >& pts,
		    const double& x, const double& y) const
{
  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      if ( pts.size() != 2 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	{
	  if ( pts[ 0 ] . size() != 2 || pts[ 1 ] . size() != 2 )
	    {
	      cerr << "Erreur sur le nombre de points pour l'interpolation"
		   << endl;
	      return x;
	    }
	  else
	    return interLin2D( pts, x, y );
	}
    case Cubique:
      if ( pts.size() != 4 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	{
       	  if ( pts[ 0 ] . size() != 4 || pts[ 1 ] . size() != 4 ||
       	       pts[ 2 ] . size() != 4 || pts[ 3 ] . size() != 4 )
       	    {
       	      cerr << "Erreur sur le nombre de points pour l'interpolation"
       		   << endl;
       	      return x;
       	    }
       	  else
       	    return interCub2D( pts, x, y );
	}
    default:
      cerr << "Fonction inconnue '" << Fonction << "'" << endl;
      return 0.;
    }
}

/**
 * Méthode permettant d'interpoler la valeur du bruit en trois
 * dimensions. Le bruit est donné par plusieurs éléments de la
 * permutation (générés par bruit).
 *
 * Pour une interpolation linéraire ou cosinus, il faut deux points
 * (l'entier précédent et l'entier suivant), pour une interpolation
 * cubique, il faut utiliser quatre points.
 *
 * @param pts valeur des bruits pour les entiers entourant la valeur
 * pour laquelle on veut la valeur du bruit.
 * @param x valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 * @param y valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 * @param z valeur pour laquelle on veut le bruit, ramenée dans
 * l'intervalle [0;1].
 * @return la valeur interpolée du bruit.
 */
double Bruit::inter(const vector< vector< vector< double > > >& pts,
		    const double& x, const double& y, const double& z) const
{
  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      if ( pts.size() != 2 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	{
	  if ( pts[ 0 ] . size() != 2 || pts[ 1 ] . size() != 2 )
	    {
	      cerr << "Erreur sur le nombre de points pour l'interpolation"
		   << endl;
	      return x;
	    }
	  else
	    {
	      if ( pts[ 0 ][ 0 ] . size() != 2 || pts[ 0 ][ 1 ] . size() != 2 ||
		   pts[ 1 ][ 0 ] . size() != 2 || pts[ 1 ][ 1 ] . size() != 2 )
		{
		  cerr << "Erreur sur le nombre de points pour l'interpolation"
		       << endl;
		  return x;
		}
	      else
		return interLin3D( pts, x, y, z );
	    }
	}
      break;
    case Cubique:
      if ( pts.size() != 4 )
	{
	  cerr << "Erreur sur le nombre de points pour l'interpolation"
	       << endl;
	  return x;
	}
      else
	{
	  if ( pts[ 0 ] . size() != 4 || pts[ 1 ] . size() != 4 ||
	       pts[ 2 ] . size() != 4 || pts[ 3 ] . size() != 4 )
	    {
	      cerr << "Erreur sur le nombre de points pour l'interpolation"
		   << endl;
	      return x;
	    }
	  else
	    {
	      if ( pts[ 0 ][ 0 ] . size() != 4 || pts[ 0 ][ 1 ] . size() != 4 ||
		   pts[ 0 ][ 2 ] . size() != 4 || pts[ 0 ][ 3 ] . size() != 4 ||
		   pts[ 1 ][ 0 ] . size() != 4 || pts[ 1 ][ 1 ] . size() != 4 ||
		   pts[ 1 ][ 2 ] . size() != 4 || pts[ 1 ][ 3 ] . size() != 4 ||
		   pts[ 2 ][ 0 ] . size() != 4 || pts[ 2 ][ 1 ] . size() != 4 ||
		   pts[ 2 ][ 2 ] . size() != 4 || pts[ 2 ][ 3 ] . size() != 4 ||
		   pts[ 3 ][ 0 ] . size() != 4 || pts[ 3 ][ 1 ] . size() != 4 ||
		   pts[ 3 ][ 2 ] . size() != 4 || pts[ 3 ][ 3 ] . size() != 4 )
		{
		  cerr << "Erreur sur le nombre de points pour l'interpolation"
		       << endl;
		  return x;
		}
	      else
		return interCub3D( pts, x, y, z );
	    }
	  break;
	}
    default:
      cerr << "Fonction inconnue '" << Fonction << "'" << endl;
      return 0.;
    }     
      
}

/**
 * Produit la valeur du bruit unidimensionelle.
 *
 * On génère la liste des points nécessaire à l'interpolation et
 * normalise l'argument par rapport à la période.
 *
 * @param x valeur de la "graine".
 *
 * @return le bruit interpolé entre les éléments de la permutation.
 */  
double Bruit::bruitInterpole(const double& x) const
{
  int i( floor( abs( x ) / dPeriode ) );

  vector< double > points;

  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      points = list_of( bruit( i ) )( bruit( i + 1 ) );
      break;
    case Cubique:
      points = list_of( bruit( i - 1 ) )( bruit( i ) )( bruit( i + 1 ) )
	( bruit( i + 2 ) );
      break;
    }

  return inter( points, modulo( x, Periode ) / dPeriode );
}

/**
 * Produit la valeur du bruit bidimensionelle.
 *
 * On génère la liste des points nécessaire à l'interpolation et
 * normalise les arguments par rapport à la période.
 *
 * @param x valeur de la "graine" 1.
 * @param y valeur de la "graine" 2.
 *
 * @return le bruit interpolé entre les éléments de la permutation.
 */  
double Bruit::bruitInterpole(const double& x, const double& y) const
{
  int i( floor( abs( x ) / dPeriode ) ), j( floor( abs( y ) / dPeriode ) );
  i /= Periode ;
  j /= Periode ;
  
  vector< vector< double > > points;

  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      points = list_of< vector < double > >( list_of( bruit( i, j ) )
					     ( bruit( i + 1, j ) ) )
	( list_of( bruit( i, j + 1 ) )( bruit( i + 1, j + 1) ) );
      break;
    case Cubique:
      points =
	list_of< vector < double > >( list_of( bruit( i - 1, j - 1 ) )
				      ( bruit( i, j - 1 ) )
				      ( bruit( i + 1, j - 1 ) )
				      ( bruit( i + 2, j - 1 ) ) )
	( list_of( bruit( i - 1, j ) )( bruit( i, j ) )
	  ( bruit( i + 1, j ) )( bruit( i + 2, j ) ) )
	( list_of( bruit( i - 1, j + 1 ) )( bruit( i, j + 1 ) )
	  ( bruit( i + 1, j + 1 ) )( bruit( i + 2, j + 1 ) ) )
	( list_of( bruit( i - 1, j + 2 ) )( bruit( i, j + 2 ) )
	  ( bruit( i + 1, j + 2 ) )( bruit( i + 2, j + 2 ) ) );
      break;
    }

  return inter( points,
		modulo( x, Periode ) / dPeriode,
		modulo( y, Periode ) / dPeriode );
}

/**
 * Produit la valeur du bruit tridimensionelle.
 *
 * On génère la liste des points nécessaire à l'interpolation et
 * normalise les arguments par rapport à la période.
 *
 * @param x valeur de la "graine" 1.
 * @param y valeur de la "graine" 2.
 * @param z valeur de la "graine" 2.
 * @return le bruit interpolé entre les éléments de la permutation.
 */  
double Bruit::bruitInterpole(const double& x, const double& y,
			     const double& z) const
{
  int i( floor( abs( x ) / dPeriode ) ), j( floor( abs( y ) /dPeriode ) ),
    k( floor( abs( z ) / dPeriode ) );
  i /= Periode ;
  j /= Periode ;
  k /= Periode;
  
  vector< vector< vector< double > > > points;
  vector< vector< double > > Plan;

  switch ( Fonction )
    {
    case Lineaire:
    case Cosinus:
      Plan = list_of< vector < double > >( list_of( bruit( i, j, k ) )
					   ( bruit( i + 1, j, k ) ) )
	( list_of( bruit( i, j + 1, k ) )( bruit( i + 1, j + 1, k ) ) );
      points.push_back( Plan );
      Plan.clear();
      Plan = list_of< vector < double > >( list_of( bruit( i, j, k + 1 ) )
					   ( bruit( i + 1, j, k + 1 ) ) )
	( list_of( bruit( i, j + 1, k + 1 ) )( bruit( i + 1, j + 1, k + 1 ) ) );
      points.push_back( Plan );     
      break;
    case Cubique:
	Plan = list_of< vector < double > >
	( list_of
	  ( bruit( i - 1, j - 1, k - 1 ) )
	  ( bruit( i    , j - 1, k - 1 ) )
	  ( bruit( i + 1, j - 1, k - 1 ) )
	  ( bruit( i + 2, j - 1, k - 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j, k - 1 ) )
	  ( bruit( i    , j, k - 1 ) )
	  ( bruit( i + 1, j, k - 1 ) )
	  ( bruit( i + 2, j, k - 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 1, k - 1 ) )
	  ( bruit( i    , j + 1, k - 1 ) )
	  ( bruit( i + 1, j + 1, k - 1 ) )
	  ( bruit( i + 2, j + 1, k - 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 2, k - 1 ) )
	  ( bruit( i    , j + 2, k - 1 ) )
	  ( bruit( i + 1, j + 2, k - 1 ) )
	  ( bruit( i + 2, j + 2, k - 1 ) )
	  );
      points.push_back( Plan );
      Plan.clear();
      Plan = list_of< vector < double > >
	( list_of
	  ( bruit( i - 1, j - 1, k ) )
	  ( bruit( i    , j - 1, k ) )
	  ( bruit( i + 1, j - 1, k ) )
	  ( bruit( i + 2, j - 1, k ) )
	  )
	( list_of
	  ( bruit( i - 1, j, k ) )
	  ( bruit( i    , j, k ) )
	  ( bruit( i + 1, j, k ) )
	  ( bruit( i + 2, j, k ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 1, k ) )
	  ( bruit( i    , j + 1, k ) )
	  ( bruit( i + 1, j + 1, k ) )
	  ( bruit( i + 2, j + 1, k ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 2, k ) )
	  ( bruit( i    , j + 2, k ) )
	  ( bruit( i + 1, j + 2, k ) )
	  ( bruit( i + 2, j + 2, k ) )
	  );
      points.push_back( Plan );
      Plan.clear();
      Plan = list_of< vector < double > >
	( list_of
	  ( bruit( i - 1, j - 1, k + 1 ) )
	  ( bruit( i    , j - 1, k + 1 ) )
	  ( bruit( i + 1, j - 1, k + 1 ) )
	  ( bruit( i + 2, j - 1, k + 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j, k + 1 ) )
	  ( bruit( i    , j, k + 1 ) )
	  ( bruit( i + 1, j, k + 1 ) )
	  ( bruit( i + 2, j, k + 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 1, k + 1 ) )
	  ( bruit( i    , j + 1, k + 1 ) )
	  ( bruit( i + 1, j + 1, k + 1 ) )
	  ( bruit( i + 2, j + 1, k + 1 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 2, k + 1 ) )
	  ( bruit( i    , j + 2, k + 1 ) )
	  ( bruit( i + 1, j + 2, k + 1 ) )
	  ( bruit( i + 2, j + 2, k + 1 ) )
	  );
      points.push_back( Plan );
      Plan.clear();
      Plan = list_of< vector < double > >
	( list_of
	  ( bruit( i - 1, j - 1, k + 2 ) )
	  ( bruit( i    , j - 1, k + 2 ) )
	  ( bruit( i + 1, j - 1, k + 2 ) )
	  ( bruit( i + 2, j - 1, k + 2 ) )
	  )
	( list_of
	  ( bruit( i - 1, j, k + 2 ) )
	  ( bruit( i    , j, k + 2 ) )
	  ( bruit( i + 1, j, k + 2 ) )
	  ( bruit( i + 2, j, k + 2 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 1, k + 2 ) )
	  ( bruit( i    , j + 1, k + 2 ) )
	  ( bruit( i + 1, j + 1, k + 2 ) )
	  ( bruit( i + 2, j + 1, k + 2 ) )
	  )
	( list_of
	  ( bruit( i - 1, j + 2, k + 2 ) )
	  ( bruit( i    , j + 2, k + 2 ) )
	  ( bruit( i + 1, j + 2, k + 2 ) )
	  ( bruit( i + 2, j + 2, k + 2 ) )
	  );
      points.push_back( Plan );
      Plan.clear();
      break;
    }
  
  return inter( points,
		modulo( x, Periode ) / dPeriode,
		modulo( y, Periode ) / dPeriode,
		modulo( z, Periode ) / dPeriode );
}

/**
 * Méthode qui réalise l'interpolation linéaire unidimentionelle entre
 * deux points.
 *
 * On fait ici l'interpolation linéaire et avec un polynome en cos(x),
 * car le nombre de points nécessaires et la méthode sont les mêmes.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @return valeur intermédiaire entre les deux points donnés.
 */
double Bruit::interLin1D(const vector< double >& pts,
			 const double& x) const
{
  if ( abs( x - 1 ) < 1.e-7 )
    return pts[ 0 ];

  switch ( Fonction )
    {
    case Lineaire:
      return pts[ 0 ] * ( 1 - x ) + pts[ 1 ] * x;
    case Cosinus:
      {
	double k( ( 1. - cos( x * Pi ) ) / 2. );
	return pts[ 0 ] * ( 1 - k ) + pts[ 1 ] * k;
      }
    case Cubique:
      cerr << "Ne devrait pas apparaître, erreur de logique : " << endl
	   << "Une interpolation cubique est demandée dans Bruit::interLin1D"
	   << endl;
      return 0.;
    default:
      cerr << "Fonction inconnue : '" << Fonction << "'" << endl;
      return 0.;
    }
}

/**
 * Méthode qui réalise l'interpolation linéaire bidimentionelle entre
 * deux points.
 *
 * On interpole en premier selon un axe, puis on interpole entre les
 * deux résultats sur le deuxième axe.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @param y paramètre pour lequel on veut la valeur intermédiaire.
 * @return la valeur intermédiaire entre les quatre points donnés.
 */
double Bruit::interLin2D(const vector< vector< double > >& pts,
			 const double& x, const double& y) const
{
  double i1( interLin1D( pts[ 0 ], x ) );
  double i2( interLin1D( pts[ 1 ], x ) );

  vector< double > tmp = list_of( i1 )( i2 );

  return interLin1D( tmp, y );
}

/**
 * Méthode qui réalise l'interpolation linéaire tridimentionelle entre
 * deux points.
 *
 * On interpole en premier sur deux plans, puis on interpole entre les
 * deux résultats sur le dernier axe.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @param y paramètre pour lequel on veut la valeur intermédiaire.
 * @param z paramètre pour lequel on veut la valeur intermédiaire.
 * @return la valeur intermédiaire entre les huit points donnés.
 */
double Bruit::interLin3D(const vector< vector< vector< double > > >& pts,
			 const double& x, const double& y, const double& z)
  const
{
  double i1( interLin2D( pts[ 0 ], x, y ) );
  double i2( interLin2D( pts[ 1 ], x, y ) );

  vector< double > tmp = list_of( i1 )( i2 );

  return interLin1D( tmp, z );
}

/**
 * Méthode qui réalise l'interpolation cubique unidimentionelle entre
 * quatre points.
 *
 * On a besoin ici de quatre points car on veut déterminer les quatres
 * coefficients de la cubique passant par ces points.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @return valeur intermédiaire entre les deux points donnés.
 * 
 */
double Bruit::interCub1D(const vector< double >& pts,
			 const double& x) const
{
  if ( abs( x - 1 ) < 1.e-7 )
    return pts[ 1 ];

  double A3( pts[ 3 ] - pts[ 2 ] + pts[ 1 ] - pts[ 0 ] );
  double A2( pts[ 0 ] - pts[ 1 ] - A3 );
  double A1( pts[ 2 ] - pts[ 0 ] );
  double A0( pts[ 1 ] );
  
  return A3 * pow( x, 3 ) + A2 * pow( x, 2 ) + A1 * x + A0;
}

/**
 * Méthode qui réalise l'interpolation cubique bidimentionelle entre
 * quatre ensembles de quatres points.
 *
 * On interpole en premier selon un axe, puis on interpole entre les
 * quatres résultats sur le deuxième axe.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @param y paramètre pour lequel on veut la valeur intermédiaire.
 * @return la valeur intermédiaire entre les seize points donnés.
 */
double Bruit::interCub2D(const vector< vector< double > >& pts,
			 const double& x, const double& y) const
{
  double i1( interCub1D( pts[0], x ) );
  double i2( interCub1D( pts[1], x ) );
  double i3( interCub1D( pts[2], x ) );
  double i4( interCub1D( pts[3], x ) );
  
  vector< double > tmp = list_of( i1 )( i2 )( i3 )( i4 );
   
  return interCub1D( tmp, y );
}

/**
 * Méthode qui réalise l'interpolation cubique tridimentionelle entre
 * soixante-quatre points.
 *
 * On interpole en premier sur quatre plans, puis on interpole entre les
 * quatre résultats sur le dernier axe.
 *
 * @param pts conteneur des points utilisés pour l'interpolation.
 * @param x paramètre pour lequel on veut la valeur intermédiaire.
 * @param y paramètre pour lequel on veut la valeur intermédiaire.
 * @param z paramètre pour lequel on veut la valeur intermédiaire.
 * @return la valeur intermédiaire entre les soixante-quatre points
 * donnés.
 */
double Bruit::interCub3D(const vector< vector< vector< double > > >& pts,
			 const double& x, const double& y, const double& z)
  const
{
  double i1( interCub2D( pts[0], x, y ) );
  double i2( interCub2D( pts[1], x, y ) );
  double i3( interCub2D( pts[2], x, y ) );
  double i4( interCub2D( pts[3], x, y ) );
   
  vector< double > tmp = list_of( i1 )( i2 )( i3 )( i4 );
   
  return interCub1D( tmp, z );
}

/**
 * Permet d'obtenir le "modulo" d'un nombre réel par un nombre entier
 * (qui est positif). L'opération revient à effectuer :
 *
 * \f[ \left( \lfloor x \rfloor ~\% n ~\right) + \left( x - \lfloor x
 * \rfloor \right) \f]
 *
 * @param nbr nomre réel que l'on veut ramener dans l'intervalle
 * voulu.
 * @param base nombre entier, limite supérieure de l'intervalle
 * construit.
 * @return la valeur qui corespond à \b \e nbr dans l'intervalle \f$
 * \left[ 0, \mathrm{base} \right] \f$.
 */
double Bruit::modulo(const double& nbr, const unsigned& base) const
{
  double n( nbr );//, ret( base );

  if ( n < -1.e-7 )
    n *= -1.;

  while ( n > base )
    n -= base;
    
  return n;
}

// /**
//  * Permet d'élever un nombre entier à une puissance entière, cette
//  * fonction n'existant pas dans la librairie standard.
//  *
//  * Pour aller plus vite, on multiplie à chaque itération par \f$
//  * \mathrm{mantisse}^2 \f$.
//  *
//  * @param mantisse nombre entier.
//  * @param exposant puissance à laquelle on veut élever \c mantisse.

//  * @return \f$ \mathrm{mantisse}^{\mathrm{exposant}} \f$.
//  */
// int Bruit::puissance(const int& mantisse, const unsigned& exposant) const
// {
//   return puissance< 2u >( mantisse );
//   /*
//   int ret( 1 );
//   if ( exposant % 2 == 1 )
//     ret = mantisse;

//   for (unsigned i(0); i < exposant / 2; i++)
//     ret *= mantisse * mantisse;

//   return mantisse;
//   */
// }
